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ABSTRACT 

A digital computer program was developed to find the buckling 
coefficient for rectangular plates with all edges simply supported or 
with all edges clamped. A finite difference technique is used to solve 
the partial differential equation for the deflection of a plate classi- 
cally treated as having only a small deflection compared to its thick- 
ness. The program was prepared to take for an input the stresses at 
nodes formed by grids dividing the plate into rectangles. The stresses 
and deflections at each node are used in writing difference equations. 
An extrapolation formula is featured in the program which allows a 
close approximation to the buckling coefficient without necessitating 
the use of a large number of grid nodes. The program was written in 
FORTRAN 60 but must be used as a FORTRAN 63 program to take advantage 
of some of its inherent flexibilities. Information is provided in 
the output of the program which aids in evaluating the reliability of 


a solution. 
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18 Introduction. 

The determination of the initial buckling load of homogeneous plates 
has been a subject of interest for many years. A very good compilation 
of past approaches may be found in the work by Gerard and Becker [8]. 

The methods used were varied and, in general, they may be classified in- 
to two categories: (1) variational methods (commonly known as energy 
methods when applied to mechanics), the best known of which are those 
attributed to Ritz and B. G. Galerkin and (2) numerical methods in which 
finite difference approximations to the partial differential equations 
of the deflection of the plate at a sufficient number of points on the 
plate result in a set of algebraic equations. The methods under the 
first category result in approximating an infinite set of infinite 
series equations by a finite set of equations which must be solved simul- 
taneously. The second method involves the coefficient matrix of a set 
of algebraic equations whose eigenvalue is found or its determinant 
evaluated for a given value of an appropriate parameter. The parameter 
is varied until the determinant vanishes. 

The first method results in a solution limited to the problem for 
which it was developed. Further, it does not promise to be simple for 
cases which may be more on the practical side. The second method was 
shown to require a considerable amount of labor and time in evaluating 
the determinant of a large matrix, even for a simple cage.” Computeriza- 
tion suggests itself as a remedy to such limitations encountered in both 

E conclusion by J. Yardley in his thesis "Applications of Finite 


Difference Equations to Buckling Problems of Rectangular Plates", 
Washington University, (1948). 


methods. 

With this in view, a digital computer program was developed to find 
the buckling coefficient for rectangular plates. The program was written 
for plates with all edges simply supported or all edges clamped. In 
writing the program the governing partial differential equation for the 
deflection of a rectangular plate based on the classical treatment of a 
thin plate having only a small deflection was used. Hence, all problems 
that are solved using the program will also be based on the same treat- 
ment. 

There are two steps that must be accomplished to get the initial 
buckling load. First, the governing partial differential equation is 
replaced by a set of linear algebraic equations by approximating it with 
difference quotients term by term. Second, the highest eigenvalue of 
the coefficient matrix of the set of equations must be found. The buck- 
ling coefficient is inversely proportional to the eigenvalue of the 


matrix. This will be shown in the following section. 


2. Mathematical Basis. 
The governing partial differential equation of the deflection of a 
plate in equilibrium in the absence of body forces and which is under 


the action of forces in its middle plane tee 


iw 4 2 Aw aw = DW му дм + 2м Su) 1 
axt ax dy” " y^ р TÉ y? Say e 


anime keno. S. P., and J. M. Gere, Theory of Elastic Stability 
(New York, Toronto, London: McGraw-Hill Book Company, 1961), p. 348. 


where 
w = lateral deflection of the plate 


x,y = cartesian coordinates. These are conveniently taken 
parallel to the sides of the plate. 


D = flexural rigidity of the plate 
N = the normal force per unit width in the x-direction 


N = the normal force per unit width in the y-direction 


На the shearing force per unit width in the x and y directions 


The forces М в, М n^ may be expressed in terms of a common parameter 


y’ 
N. Thus we put 


= 


„Troy 


N, = N Е (х,у) 
n. -N n (қыт) 

where Ffy and F. are specified functions of x and y. Assuming that 
such a relationship exists, the problem of finding the buckling coef- 
ficient becomes a problem of finding the smallest value of N that would 
cause the plate to start to buckle, or the critical value of N. To do 
this using finite difference technique, Eqn. 1 must be replaced first 
with a set of linear algebraic equations. 

The plate is divided into integral parts in the x and y directions 
so that rectangular meshes of uniform size are formed. This manner of 
dividing the plate is advantageous because it is not limited by the 
aspect ratio of the plate and thus makes programming easier. At each 
corner, or "node", of the meshes a finite difference equation is writ- 
ten in approximation to Eqn. 1. The set of algebraic equations formed 
may be written in matrix form as 


Абы = kBw (2) 


where 
pe = square coefficient matrix for the left side of Eqn. l. 
B = square coefficient matrix for the right side of Eqn. l. 
w = column vector of the plate lateral deflection at nodes. 


k 


function of the parameter N. 


Eqn. 2 may be rewritten as 


(c -inw-o (3) 


where 


I = the identity matrix of the same order as C. 
Except for the trivial case when w is zero, Eqn. 3 is true only when 
x= ). 
where А is an eigenvalue of the matrix C. Since we are seeking the 
smallest value of k, the highest eigenvalue A must be found. The de- 
tailed mathematical formulation for the matrices is discussed in Ap- 
pendix I. 

It was realized that theoretically the outlined procedure of find- 
ing the buckling coefficient will give a value close to the correct buck- 
ling coefficient only when the number of points used is large. There 
are a number of objections to satisfying such a requirement. Dealing 
with a large number of points will require a large amount of computer 
storage since for a given square matrix a storage equal to the square of 
its order is required. The second objection to the use of a large number 
of pofi#its in the plate is that too much computer time will be involved. 
The last objection is that there is possibility of round-off which would 


introduce an error unless computation using double precision arithmetic 


10 


is resorted to. 

To eliminate the necessity of using a large number of points on the 
plate an extrapolation formula is used. In doing this one effectively 
looks for an estimate of the error inherent in the finite difference 
approximations of the governing partial differential equation when using 
a finite number of points in the plate. In other words the inherent 
error is the result of replacing an infin ites ina} quantity with a finite 


one. For this purpose the error € was assumed to take the form 
2 4 2. 2 
€ = K, H, tKa Hx D K, Hy " Ks Hy" v Ke ж Hy 


Extrapolation was achieved by solving simultaneously for К, Six equa- 


tions of the form 


К = К+ К.Н + Kale +Ky Hy + Ks Hy + к. ННЯ 


where 


K = the computed buckling coefficient for a given choice of 
mesh 


K, = the extrapolated value of the buckling coefficient 


K, ,K,sK, K; Kç = constants 


H. = the length of the mesh used 
е. = the width of the mesh 
A more detailed discussion of the estimate of the error and extra- 


polation for the buckling coefficient is given in Appendix I. 


3. Cases Considered. 

Several problems in plate buckling with known solutions were solved 
using the program in order to compare the results that were obtained 
with it. The discussion below will indicate the extent of agreement with 


11 


previous solutions. All cases will be discussed with reference to tables 
which are the final output of the program. Tables I-VII may be found at 
the end of this section on pages 21-27., The information provided by the 
tables is 
Aspect Ratio - length/width 

K, = extrapolated value of the buckling coefficient. 

K,,K,,K,,K.,K. - constants calculated for correcting K. 

K = the computed value of the buckling coefficient 
= corrected value of K. 


Moros 
K, and Қайы. c caq should agree very closely if round-off errors are 
not excessive since they are both based on the same equations. All 
quantities are dimensionless. 
The analytical solutions to the problems that were considered were 
cast in the form 
N K TT2p/b2 
where 
b = width of the plate 
K' = a dimensionless constant 
In the program K was computed when the above formula is written in the 
form 
N = KD/b2 


Comparison of solutions wíll be made between K, and к'тї^. 


1 


CASE I. Simply Supported Plate Under Uniform Compression in the X- 


Direction. 


For this case the common parameter is N - N, so that 


12 


F =F = 0 
2) ху 


The formula for the critical load of rectangular plates with aspect 


ratio = a/b where a is the length of the plate and b is its width 2 


2. 2. 
N° TD (2 + b.) 
a b> b О, 

Еог а/Ъ = 3/4 we have 

N_.. = 42.836826 D/b” 
Entering Table I we get 

K, = 42.836819 
Comparing this with the coefficient of the theoretical critical load 
we get a 

Difference = .000016% 
This small difference substantiates the form of the extrapolation formula 
assumed. 

That the finite difference solution is accurate in this case may be 
shown also by considering the eigenvectors using the set of points when 
there are 10 divisions in the x-direction and 12 divisions in the y- 
direction. It is known that the deflection surface of the buckled plate 


may be represented by the question: 


wai, sin(Wx/a) sin({y¥/b) (4) 
where 
841 = constant 
Referring to the figure in Table I and setting the deflection w = l at 
point #1, aj, can be determined. We get 
ауу = 12.503206 
9 
Ibid., p.353. 
4 


Ibid., p. 327. 13 


With 241 known we get for point #46 


WAG ~ 12.503206 


The computer solution for the eigenvectors gives for point #46 


W, 6? computer - 12.502909 
Comparing the values obtained we have 


Difference = .0023% 


CASE II. Simply Supported Plate Subjected to Pure Bending. 
The figure in Table II shows the plate under pure bending load. 


The parameter chosen is the maximum stress designated as N in the figure 


so that 


Hj 
II 


1 - 2y/b 
F = O 
Е =0 
The solution to this problem using three equations of an infinite set 
is given by Timoshenko Bec 
No = 24.1117 D/b* 
= 239 D/b? 
Entering Table II we get for an extrapolated value of the buckling co- 
efficient 
K, = 238 


Comparing this with the constant for the critical load we get a 


Difference = .4% г 


?Ibid., p.377. 
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CASE III. Simply Supported Plate Under Compressive Load Varying 
Linearly in the Direction of Loading. 


Referring to the figure in Table III and choosing N as the common 


parameter we have 


F = 2x/ae 
x 
F =O 
y 
E. = 2y/a - b/a 


The solution to this by Libove et. "E using matrix iteration methods 
on the matrix obtained from replacing an infinite set of equations is 

н _ = 2.92 TV2p/b2 

сї 
- 28.8 D/b? 

Going to Table III we get for a coefficient 

Ki z 27.9 
À comparison between the two coefficients gives a 

Difference = 3% 
This difference is very much greater than for the other cases that have 


been considered. This may be due to the presence of shear stress whose 


effects are discussed in the next case. 


CASE IV. Simply Supported Plate Subjected to Pure Shear. 


In Table IV is shown a plate with aspect ratio a/b = 2.5 subjected 
to pure shear. For this case the common parameter is N = Ұқ which gives 


the following: 


6 

Libove, C., S. Ferdman, and J. J. Reusch, "Elastic Buckling of a 
Simply Supported Plate Under a Compressive Stress that Varies Linearly 
in the Direction of Loading", NACA TN no. 1891, p. 18, (1949). 
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F =1 
xy 
Timoshenko and Gere obtained the sorteo to this problem by replacing 
an infinite set of equations with five equations and equating to zero 


- 


their determinant. This gave a result of 


н. = 6.112? 
= 60 рљ? 

Going into Table IV we get a coefficient 

K; = 60 
On the basis of the available significant figures in the published value, 
the solutions have a 

Difference = 0% 
it will be noted that only three computed values are presented in Table 
IV. Originally six computed values of the buckling coefficient were 
used. These values indicated that the true buckling coefficient had a 
value lower than that obtained with the grid choice that gave the greatest 
number of interior nodes. The extrapolated coefficient had a value that 
was practically twice the value of any of the six results - which made 
it hard to accept as correct. An examination of the eigenvalue sub- 
routine output revealed that the traces of the C-matrices for two grid 
choices were non-zero. This is contrary to what is expected for this 
case according to theory which will not be discussed here. For another 
grid choice it was found that the iteration for the eigenvalue had not 


converged. The three computed buckling coefficients based on these grid 


choices were all deemed unreliable. As a remedy to the situation, the 


/rimoshenko, S. P. and J. M. Gere, loc. cit., p. 382. 
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extrapolation formula was used with only the second érder comes. - 
This required only three separate results, which were available. It 
will be noted that the extrapolation gives adequate accuracy for 2- 

significant-figure-comparison. 

In Case III where shear is combined with compressive stress the 
effect of shear is made evident in the extrapolated value of the buck- 
ling coefficient which is closest to the computed value using the coars- 
est mesh. The traces and eigenvalues do not give enough information to 
use as a basis for discarding any data as unreliable. However, discard- 
ing the result for the coarsest mesh and extrapolating without the 


2 


A, Hy term gives a value K, = 28.6 which compares with the published 


1 


value to an agreement of .7%. 


CASE V. Clamped Plate Under Uniform Compression. 

This is the same problem treated in Case I with the edge condition 
changed. All the stresses except N. will be zero. S. Levy gave a solu- 
tion to this problem based on an asymptotic approximation for an in- 
finite determinant." He cited other values that compared to his solu- 
tion by 2-975. He gives 

N = 11.659 15, 7? 
= 115.07 D/b* 
The computer solution to this problem, from Table V, gives a coefficient 
of 
к, = 115.41 


Comparing solutions we have a 


Difference = .3% 


Say, S., "Buckling of Rectangular Plates With Built-In Edges", 
Journal of Applied Mechanics, Vol. 9, pp. А171-А174, (1942). 
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CASE VI. Square Plate With Clamped Edges Sub jected to Hydrostatic 


Compression and Constant Shear. 


Taking the shear force as the common parameter we have 


F =F = 1.5 


The published value of the buckling coefficient w. 
N _ = 3.244Y(2p/b2 
cr 
= 32.0 р/ь? 
Table VI gives the computer solution to this problem as having a co- 
efficient of 
K; = 32.3 


À comparison between coefficients gives a 


Difference = .97 


CASE VII. Clamped Plate Subjected to Pure Shear. 


In this problem all the stresses are zero except shear so that 


F =F =O 
x y 


F ç =1 
xy 
and the common parameter is the shear force itself.  Interpolating from 
curves drawn by B. Budiansky and R. Connor in their solution to the same 
problem we Obras lo 
N _ = 9.9%(2p/b2 
cr 
- 97 р/ь? 


anos S. P. and J. M. Gere, loc. cit., p. 386. 


10; udiansky, B. and R. W. Connor, "Buckling Stresses of Clamped 
Rectangular Flat Plates in Shear". NACA TN No. 1559, p. 10, (1948). 
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Going into Table VII we obtain for a computer solution a coefficient of 
z 96 
K, 9 
À comparison between the computer solution and the value obtained from 


the curves gives a 


Difference = 1% , 


4. Conclusions and Recommendations. 

This program is capable of solving a variety of problems in buck- 
ling of rectangular plates that have simply supported or clamped edges. 
In the only case where an exact solution was available the extrapolated 
value of the buckling coefficient agreed with it to within .000016%. In 
other cases considered it was demonstrated that the solutions obtained 
using the program compare very closely to published solutions. 

The information provided by the program output is considered ade- 
quate. The trace of the C-matrix was found useful, as exemplified in 
the cases solved involving pure shear. In any problem, the computed 
buckling coefficients give the user a rough idea of what the magnitude 
of the extrapolated buckling coefficient K, should be. This was put to 
use in all the cases considered where shear stress was involved. It will 
be profitable to evaluate the reliability of a result by studying the 
iterates and number of iterations which are part of the eigenvalue sub- 
routine output. If there are 16 iterations there is a possibility that 
the iteration has not converged. In this event the user should examine 
the iterates closely. 

Of the many edge conditions involved in problems of plate buckling, 
only two were considered. The program will prove more useful if this 


feature is extended to other combinations. 
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It is felt that the treatment of problems involving shear is in- 
adequate. A square plate subjected to hydrostatic pressure was divided 
into 6 divisions in the x-direction and 12 divisions in the y-direction. 
It was observed that the trace of the C-matrix for such a grid choice 
agreed to the 9th significant figure with the trace of the C-matrix when 
the order of dividing the plate was reversed, i.e., 12 divisions in the 
x-direction and 6 divisions in the y-direction. 

For a square plate subjected to hydrostatic pressure and shear the 
agreement is only up to the 3rd significant figure. It is evident that 
such a behaviour is due to shear, but the explanation has not been found. 

The method of approximating the partial derivatives on the right 
side of the governing equation effectively required the use of a mesh 
twice the size of the selected mesh as far as the mixed partial deriva- 
tive term is concerned. Since it is this term that is associated with 
shear, it might be one reason for the strange behaviour of the program 
when shear is involved. It is possible that the method suggested in Ap- 
pendix I of using Green's Theorem to approximate the right side of the 
governing partial differential equation will remedy the difficulty en- 
countered with shear. This method also makes it possible to make the 
C-matrix symmetric and, as a consequence, the problem of slow convergence 
may also be solved, since there are many eigenvalue subroutines which 


can handle symmetric matrices effectively. 
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APPENDIX I 


MATHEMATICAL DEVELOPMENTS 


I- Finite Difference Approximations. 
The partial differential equation to be approximate is 


4 2 2 2 
= [М QW 2 N N 
wage) s 


^ 4 
v^w - QW +2 aw + OW 
D xt Әх“ду ду 
І-А. Simply Supported Case. 
The boundary conditions at all edges of the plate are 
w=0 (1-2) 
2 
V w^0 (1-3) 
where we define у? as the Laplacian operator which is 
2 2 
v^ m acie GL за 
2 2 
OR Әу in two dimensions. 
Since it is known that the values of w = 0 at the boundaries by virtue 


of B.C. (1-2) and, since we need only as many equations as there are un- 


knowng,only equations for interior points need be written. 





Fig. I-1 
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The left side of Eqn. 1-1 may be approximated easily once ^w is 
obtained. This is done with reference to Fig. I-l. Let the numbered 
points be the mesh nodes when the plate is divided into meshes with the 
mesh sides of length H. and H At point #0 using the conventional 


method of approximating the partial derivative we have 








Danes Wa + Wz, - 210 (1-4) 
9x? He 

2 
ы 2 W + wa 2Wo (I-5) 
o y Hy 

2 


V Wo =e tiia + (veajer 2 өз) (1-6) 


Àn examination of Eqn. I-6 shows that for a given point on the plate, 
four other points will be involved in writing the equation for yw. 
The deflection we at each point will have the following coefficients 


when 1/ HH) remains factored out: 


а, = a = ME 


1 8 Hy 
= = H 
a, = a, he 
8977847857 84,7 8, (1-7) 


When one of the points involved (other than point £0) lies on a 
boundary, a slight modification in obtaining the coefficients must be 
performed. If in Fig. I-1 point #4 lies оп the boundary, its contribu- 
tion in Eqn. I-7 for the point #0 is still included, but а, will not ap- 


pear in the A-matrix and w 


4 is not used in assembling the deflection 


vector. We have seen that 


2 


Hx Hy 
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In general, for n mesh points, Eqn. 1-8 may be written as 


Y = Aw (1-9) 


where 
A = ann x n symmetric matrix 


w = a column vector of lateral deflections at interior 
nodes of the plate. 


For a square mesh the deflection coefficients will be 


| а, = 1 1 = 1,2,3,& 


a = -4 

o 
Matrix A is easily assembled if we use the five-point cross for the co- 
efficients as illustrated in Fig. I-2. The number on each node represents 


the coefficient of w at that node. 


Fig. I-2 
Having obtained y? м ме аге now ready to approximate’ E To 


get Xs we simply operate on Eqn. I-9 so that 
4 2/_2 2. 
V W = ү (V w = A (Aw) LAW (1-10) 


It should be noted that vw was formulated under the condition that w = 0 
at the boundaries and vw was formulated under the condition that vw 


= Q at the boundaries also. 
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It has been shown that only vw need be developed to obtain V 


The scheme used to write down the matrix for vu will be discussed next. 


| St. ради 


2, hd. par 


\ ' і 
= rd "Paar 





Fig. I-3 


1 


Let the mesh nodes be numbered in the manner shown in Fig. ТО И 


where two adjacent horizontal mesh lines are taken as a pair and each node 
on the lines is numbered alternating between the two lines. If we write 
down the approximations to Eqn. I-4 for each point on the plate in the 
sequence that they are numbered in Fig. I-3, using a square mesh a co- 
efficient matrix such as shown in Fig. 1-7 results. This is the matrix 
that must be multiplied by itself to get the coefficient matrix Ae for 
wae If we write the equation for a given point i the deflection coef- 
ficient will lie in the ith row and in the ith column. The other points 
involved in the equation for point i will have a deflection coefficient 
lying i the ith row and in columns bearing their respective numbers. 


Thus it is easy to locate the deflection coefficient of any point. In 


Ume scheme is based on the work of Griffin, D. S. and R. S. Varga 


in their paper "Numerical Solution of Plane Elasticity Problems", J. 
Industrial and Applied Math., Vol. 11, pp. 1046-1061, (1963). 
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Fig. I-7 which was assembled using a square mesh, the deflection at point 
4 will have a coefficient of -4 in row 4, column 4. Points 2,3,6, and 15 
will have deflection coefficients equal to 1 in the same row and lying in 
columns 2,3,6, and 15, respectively. 

Squaring matrix A we will get the matrix for vw which is illustrat- 
ed in Fig. I-8. It is now simple to check this matrix. If we follow the 
same procedure of locating the elements in any row it can be verified 
that the matrix А? which was assembled for a square mesh agrees with the 
matrix formed using the 13-point star to approximate Ju for a square 
mesh. This star is illustrated in Fig. I-4 which shows the deflection 


coefficient for each point. 





Fig. I-4 


This completes the approximation to the left side of the governing partial 


differential equation. 


Approximations of the Right Side of Eqn. І-1. 





The first attempt to approximate the right side of Eqn. I-l was to 


apply Green's Theorem to the integrals 
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3 x Du dxdy 4 СУ = -— i зды | dy) 


gx? 


Applying the theorem will result in a symmetric coefficient matrix B. 
Since À and as are symmetric, we may then deal only with symmetric 


matrices. The C-matrix can be made symmetric also by the following steps. 


A^w = kBw 
A(Aw) = КВА (Aw) 
ie suada 
w' = KA BA Я 
- = kcw' 
where 
w' = Aw 
and Ce BA "which ís symmetric Ui [REIR 


symmetric. A symmetric C-matrix is advantageous because less storage is 
required by the digital computer program and there is a wide choice of 
dependable subroutines for finding the eígenvalues of a symmetric matrix. 
However, a symmetric B-matrix can be obtained only with the use of stresses 
at intermediate points between nodes. The resulting matrix will be more 
intricate and programming more elaborate since a problem in "bookkeeping" 
arises from the use of stresses at points other than at the point under 
consideration. For these reasons a direct approximation to each term in 
the right side of Eqn. I-l was taken using the classical method of ap- 


proximating the slope at a given point. 
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í $ 
Se T 


Fig. I-5 
Using Fig. I-5 and remembering that H. and a are constant throughout the 
entire plate, the right side of Eqn. 1-7 may be approximated term by term 


in the following manner: 


D Ww SN ( Wa=wWe _ wo Mz, | 
Dx? Hx Hx Hx 
Nx Ow = Nx (94 rW,- 2 wo | (I-11) 
LS H: 
similarly 
Ny JW = Ny (w * Ms - 2 wo) (1-12) 
dy” Hy” 


For the mixed partial derivative 





Nyy БҰТЫ Nx (we + Wg = W7 47 (1-13) 
Oxdy АНН 


For the purpose of programming, 1/(H H.) was factored out as in the de- 
velopment for A. Combining Eqns. I-11, I-12 and I-13 and remembering that 


N,N, andN are expressed in terms of N and F , F , and F we obtain 
x' y xy x” y ES 





zum all (5 (Wa +W2- 2Wo) Hy + 5 v PW; - “© He + rW WI E 
DH, Hy Hx Ну 2 


The coefficient of w at each node involved may be summarized into 


oe -2(ғ н//н, + F H/H) 
DEUS F H/H, 
DE = FAM 
bs, =b, = - Fey”? 


b, =b_ =F /2 
xy 


It has been shown that the right side of Eqn. I-1 may be written in 


the form 


Right side - aid bw, қом, мі b. Wa + beWe 
DH, Hy 


In general, for n interior mesh nodes we will have the matrix equation 


Right side = kBw 
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where 
В = ап n x n coefficient matrix 
k = N/D 
An example of matrix B is shown in Fig. I-9. This corresponds to the 


plate illustrated in Fig. I-8. The loading is 


F =F = -1 
x y 
F = 2 
xy 
The common parameter in this case N = N, = - мш This completes the 


approximation to the right side of Eqn I-l. 

It will be noted that the pairing of mesh lines in Fig. I-3 depends 
on the divisions only in the y-direction. It is readily seen that the 
matrices will be different for an odd number than for an even number of 
divisions along the y-direction. The matrices for a plate divided into an 
even number of divisions along the y-direction have been illustrated in 
Figs. I-7, 1-8, and 1-9. Figs. I-10, I-11, and I-12 illustrate the 


respective matrices for an odd number of divisions in the y-direction. 


Symmetrical Cases 


There are cases where the load and deflection are symmetrical with 
respect to the geometric axes of symmetry of the plate. When this is 
known (or shown by eigenvectors obtained using the whole plate) we may 
choose to use only one quadrant of the plate. This allows use of a 
finer mesh without a corresponding increase in computer storage require- 
ments. It is the intention here that only the case when the midpoint of 
the plate has the maximum deflection will be handled by the program.  How- 


ever, where the buckled surface assumes a complete cycle of a sine curve, 


37 





say in the x-direction, then one-half of the plate may be considered and 
if it satisfies the conditions of symmetry as described previously, a 
quarter of this may be chosen for the program to handle. 

The procedure for assembling the matrices for a quadrant of a plate 
is the same as in the preceding section with the following modification. 
The whole plate is always divided so that the horizontal and vertical 
axes of symmetry are taken as mesh lines. Going back to Fig. I-5, suppose 
that points 0,1, and 3 are lying on the line of symmetry and points 2,5, 
and 6 are in the quadrant being analyzed. The coefficients aj» a, and a) 
are computed as before but а,, ас, апа ас аге now doubled since the follow- 


ing conditions obtain: 


as = ag اس ر‎ 
ao A ofa aui 
НО, MONS M 


Points 8, 4, and 7 are not used in assembling the matrix. Squaring A 
does not give a symmetric matrix. To make A? symmetric it has to be 
premultiplied by a diagonal matrix whose elements are proportíonal to 
the areas of the mesh regions associated with the corresponding matrix 
rows. This matrix may be designated as a,” and is used for obtaining 
matrix C for the quarter plate. 

The procedure described also holds equally well for the right hand 
side of the governing equation. B, is obtained after multiplying matrix 


1 
B generated for the quadrant of the plate with the diagonal matrix used 


in obtaining а 
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I-B. Clamped Edges. 


The boundary conditions at all edges are 


м = 0 (1-14) 
Pl 0 (1-15) 
ох 
ow = 0 (1-16) 
OM 


Approximation of the Left Side of Eqn. 1-1. 


Returning to the developments in the case of simply supported edges, 
it can be demonstrated that the development of A still holds. Since the 
first boundary condition for both cases is the same and since it is the 
only boundary condition used to assemble the coefficient matrix for vw, 
the matrices will be identical for both edge conditions. In assembling 
the coefficient matrix for Yw it will be necessary to take into account 


that du is nonzero on the boundaries. In Fig. I-6 we have 





Vw. - ia cv e (veles (e Wo) Ha (1-17) 
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Applying the conditions that at the right edge 


Wo - Vi - Wa = 0 
and 
w, = V, by virtue of B.C. I-15 
2 
Mig, 2T A wh = Maza 1-46 


When we take Vu by operating on Eqn. 1-7 there will be a term in Yu, 
which is nonzero (unlike the case of the simply Supported plate). The 
coefficient of the term will be the correction necessary to convert the 
matrix for simply supported plate into that of a clamped plate. This 


correction will be of the form 


Beorrection ires +, 
- 2m" 


Excluding the "corner" points, the correction above is good for all 
points adjacent to the right and left edges of the plate. For points 


adjacent to the top and bottom edges of the plate the correctíon will be 


4 
асоггесёіоп | 2/8, 


For the "corner" points there will be a contribution from points 
lying on the boundaries in both directions so that the correction term 


will be a combination of the corrections above or 


= & 
“correction > + 1/8, ) 


In terms of the matrices 
2 2 
A’, clamped = A’, simply supported + As 
where 


A =nxnd ; | 
Слан n diagonal matrix of correction terms 
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Fig. I-13 illustrates the coefficient matrix for Re using even numbers 
of divisions along the width of a clamped plate. 
Approximation of the Right Side of Eqn. I-l. 

It has been noted that the simply supported plate and the clamped 
plate had zero deflections along the boundaries. This leads to the con- 
clusion that the assembly of matrix B for a clamped plate will be exactly 


the same as for the simply supported plate. 


Estimate of the Error. 

The finite difference method of solving a partial differential equa- 
tion carries with it an inherent error resulting from the replacement of 
an infinitesimal quantity with one that is finite. To get a satisfactory 
answer the use of a considerable number of interior points will be neces- 
sary. When this is not a possibility one may use some form of extra- 
polation formula to get a good approximation to the real solution. 

There were two schemes tried to evaluate the error resulting from 
the finite difference approximations to Eqn. I-l. The first was to 
assume that the error took the form 

Б бум 
C = C, Hy 4% C2 Hy 

where c = the error 

CC, n. Bs = constants 
The calculated and extrapolated values of the buckling coefficient then 
have the relationship 

Ka = К, С" + Coby 4 
where 


ж 
lI 


computed buckling coefficient 


7 
i 


the extrapolated buckling coefficient 


4] 


For five computed values of K. the extrapolated value K, of the buckling 


1 
coefficient can be found. This has been done and, while the results were 
not very far from those obtained analytically, the method which will be 
described below was found to give better results. 

It has been m that the error in approximating Vw with finite 
differences using a square mesh can be expressed as 


е - K,H*+ Къ ЦА y Ka H° + mans rc e nn 
where 


en 


K = constants 
H = length and width of the mesh used. 

It has been a also that the same form of expression holds true for 
yw. The difference between the two will lie only in the constants. 
While such an expression applies only to Yu and vw it is reasonable 

to expect that the error for the governing partial differential equa- 
tion will vary with even powers of the mesh sides because of its form 
which is practically a combination of ү» апа uw. It will be assumed 


that the error will vary according to: 
2 4 2. 4 ага 
€ = K.H, + Ks Hx + Ky Hy + Ks Hy + Ke Hy Hy 
The relationship between the computed buckling coefficient and the extra- 
polated value will again take a form similar to Eqn. I-19. For six 


computed values of K. the extrapolated value K. of the buckling coef- 


1 


ficient can be determined. 


ta aceon tenn L. V. and V. I. Krylov, Approximate Methods of 
Higher Analysis (New York: Interscience Publishers, 1958), p. 196. 


scarborough, J. B., Numerical Mathematical Analysis, (Baltimore: 
The Johns Hopkins Press, 1962), p. 399. 
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APPENDIX II 


DESCRIPTION OF THE PROGRAM 


II-A. The Main Program. 

The main program follows closely the developments in Appendix I. 
Because of the regular pattern of the matrices resulting from the scheme 
of numbering the mesh nodes, the program was developed for variable orders 
of matrices whose forms depend on the number of interior mesh nodes used 
and the number of divisions in the y-direction (whether they are odd or 
even). The position of the elements along lines parallel to the main 
diagonal followed simple arithmetic regularity so that the generation of 
the elements was in terms of the diagonals they belonged to. 

Matrix A was shown to be symmetrical. This made it possible to 
generate the elements diagonally opposite across the main diagonal at 
the same time. This is not true in the case of matrix B because of the 
method in approximating it. 

The program follows the steps enumerated below. It is summarized in 
a general flow chart of the program on page 52. 

l. Generate А 

2. Square A 

3; Correot-A^ af the plate is clamped . 

4. Correct 2 if lines of symmetry are used 
5. Invert а 

6. Generate matrix B 

7. Correct B if lines of symmetry are used 


8. Calculate C - А7? 
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9. Find the eigenvalue of C 
10. Compute the buckling coefficient K 
11. Repeat all of the above to get enough K's for 
extrapolation 
12. Extrapolate 


13. Get eigenvectors if desired 


II-B. Subroutines 

Subroutine MATALG - This is available as a mathematical subroutine 
at the USNPGS Computer Facility. It has two options that suit the re- 
quirements of the main program. The first option is needed to invert 
the matrix a, The second option is used to solve the simultaneous 
equations to get the eigenvectors if desired, and to solve the correction 


equations to obtain KK K,,K,,K Ke» and K This subroutine 


ZAG corrected' 
is capable of providing an inverted matrix accurate to at least 9 signi- 
ficant decimal digits. 

Subroutine EIG3- This subroutine is also available as a mathematical 
subroutine. It is used to evaluate the eigenvalues for the matrix C. 
Since it does not solve for the eigenvectors, subroutine MATALG is called 


to provide them. The mathematical methods applied in this subroutine are 


discussed extensively in Ref. 6 . 
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APPENDIX III 


INSTRUCTIONS ON THE USAGE OF THE PROGRAM 


General. 

The program is written in Fortran 60 (F-60). However, it is used 
as an F-63 program for three reasons. First, extra storage is obtained 
with the use of control cards "RELOCOM" and "EXECUTER". Second, the 
variable dimensioning feature is necessary to get enough computed values 
of the buckling coefficient for different mesh sizes for the purpose of 
extrapolation. Third, the time for compilation and execution is reduced 
with the use of binary decks. It may be noted that one minute is saved 


when a listing for the program is not called for. 


III-A. Purpose. 

Program Buckle was written to compute the initial buckling load of 
rectangular plates and, if desired, to find the relative deflections of 
point on the plate at the start of buckling. This program is limited to 
plates with edges that are simply supported or to plates with clamped 


edges. 


III-B. Input Requirements. 

It is assumed that the stresses at every mesh node are known and 
stored as one dimensional arrays in XFORCE, FY, and FXY (these corres- 
pond to Р, Fy and нет? respectively). Іп the sample program, the 
computation of the stresses was incorporated as part of the main pro- 
gram. This part can be removed completely and introduced as a subprogram. 


For one particular manner of dividing the plate one data card is required. 
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The data card is divided into 7 fields of four characters each. These 
are reserved for the following parameters: 

1. Q = symmetry parameter. For the whole plate it is entered 
as 2.0. For a quarter of the plate it is entered as 1.0. 

2. Clamp = support parameter. For a simply supported plate 
it is entered as 1.0. For clamped plates it is entered as 2.0. 

3. MCC = the number of divisions in the y-direction. It must 
be right-justified and in fixed point. This may not be less than 6. 

4. NNR = the number of divisions in the x-direction. It must 
be right-justified and in fixed point. This may not be less than 6. 


5. AS 


the length of the plate. It must be in floating point. 


6. BS = the width of the plate. It must be in floating point. 


H" 


7. VECTOR = eigenvector option parameter. This is entered as 
1.0 when eigenvectors are needed, otherwise any other positive floating 
point number is entered. 
Since six computed values of K are needed to extrapolate for the 
buckling coefficient, six data cards must be prepared. Because the pro- 
gram is limited to handle a maximum of 99 internal mesh nodes when the 


whole plate is used, the product.(MCC - 1) (NNR - 1) may not exceed 99. 


III-C. Output of the Program. 
The program will have the following output. 
1. The trace of the matrix C. 
2, The iterates in finding the eigenvalue and the number of 
iterations. 
3. The first, second, and third derivatives of the poly- 
nomial used to approximate the determinant for a given eigenvalue. 
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All of the above are provided by subroutine EIG 3. 
4. The extrapolated value of the buckling coefficient K). - 


5. The constants used to evaluate the error, K K, Ko > K 


29 2 


and K,. 

6. The number of divisions in the x and y directions. 

7. The computed value of the buckling coefficient K for each 
choice of grid. 

8. The corrected values in (7), к ened 

The values in (4) and (8) should agree very closely if the extra- 

polation is correct. 

9. The eigenvectors are printed out when the eigenvector option 
is selected. The components are printed out starting with v, since W = 


1 in the program. The printout reads from left to right at the start 


of every line. 


III-D. Cautions to Users. 

The use of a quarter of a plate when the buckling is symmetrical has 
not been tested satisfactorily. While the matrices generated using this 
feature of the program were found correct, finding the eigenvalues for 
the matrix C required excessive computer times. Unless further tests are 
made, it is suggested that using the whole plate is a more reliable feature 
of the program to use. 

When solving problems with this program, it will be important to 
check the iterates and the number of iterations to find the eigenvalue 
for every C-matrix. The subroutine for finding the eigenvalue was writ- 
ten to accept as an eigenvalue the 16th iterate when convergence is slow. 
This means that an eigenvalue which requires 16 iterations to find is not 


55 


the true eigenvalue although it may be close to it. The computed buckling 
coefficient based on this should be considered unreliable. 

There is a theoretical basis which will not be discussed here for 
one to expect that the C-matrix for cases involving pure shear has a 
trace equal to zero. Hence, when one is dealing with such a problem 
he should be wary of data —- is based on a C-matrix with a trace that 
is nonzero. For this purpose a trace whose absolute value is equal to 
Р ТАЯ 

When the above difficulties are encountered, one may decide to 
change grids and replace the unreliable data. It may be more expedient 
to use the extrapolation formula with the fourth order terms omitted if 
enough data is left after discarding the unreliable results. In this 


event, at least three results must be available to use the extrapolation 


formula including only the second order terms. 


IV-E. Time Requirements. 
The program requires an average time of 3 minutes to compile and 
12 minutes to execute for a total of 15 minutes running time to solve a 


problem using six different grids. 
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APPENDIX IV 


SAMPLE PROGRAM 
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PROGRAM BUCKLE, DECK ASSEMBLY 


CONTROL CARDS 


PROGRAM BUCKLE 


END 


SUBROUTINE MATALG 


END 


SUBROUTINE EIG3 


END 
FINIS 
EXECUTER. 


DATA CARDS 
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A,B 
AS, BS 
AA 

BB 
CLAMP 
FNCR 
FORCE1 
FRY 
FY 


Kl 


NM 
Q 
VECTOR 


XFORCE 


LIST OF SYMBOLS 
arrays containing the elements of the coefficient matrices 
length and width of the plate, respectively 
Hy Hy 
H,/Hy 
edge condition parameter 
calculated buckling coefficient | 
corrected buckling coefficient 
array for shear stress 
array for stress in the y-direction 


number of pairs of horizontal mesh lines. (A single last 
line counts as a "pair'.) 


number of interior mesh nodes, also order of the matrices 
symmetry condition parameter 
option parameter for eigenvectors 


array for stress in the x-direction 
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25 
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20 
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221 


23 


22 


224 


225 


251 
253 


227 
226 
800 
26 
27 


280 


28 
282 
900 


Bl J-29J)=AA 

GO TO 172 
B(J-29J)=2.*AA 
CONTINUE 

DO 252 L= 19K6 
J=(L-1)*NT+4 
IF(Q-1.)250,25,250 
B(J=29J)=AA 

GO TO 252 
B(íJ-29J) 24. *AA 
CONTINUE 

GO TO 122 

00119 (-1%К1 

DO 20 К-5»МТ 
J=(L-1)*NT+K 
B(J-29+J)=AA 

DO 19 Kz3,NT 
J=(L-1)*NT+K 
В(.)».)-2)-АА 
CONTINUE 

GO TO 221 

DO 123 L= 1»K6 

DO 124 K= SoeNT 
Jz(L-1)*NT*K 

B( J—-29 J) AA 

DO 125 K=3eNT 
J=(L-1)*NT+K 

Bl JeJ—-2)=AA 
CONTINUE 

GO TO 22 

DO 23 J=1»K7»2 

B( J* 19NT*J) BB 
B(NT*J9»J*1) BB 

GO TO 26 

DO 224 Jz15K1052 
BCJ*19NT*J) BB 
BCJ*NT9»? J* 1) x BB 

DO 225 J= 1»NR 
B(K9*2*J) .,K2*J)- BB 
IF (Q-1e.)251»227»251 
DO 253 J= 1»NR 
B(K2*J9K9*2*J) BB 
GO TO 800 

DO 226 Jz1,NR 
B(K2*J)»K9*2* J.) 22, *BB 
GO TO 900 

DO 27 J= l»K7»2 
B(J+1,J)=BB 

00282 J= K5»NN»2 
ІҒ(С-1.)280»28»280 
Bt J*1»J) BB 

GO TO 282 

B(J+lsJ )=2.*BB 
CONTINUE 

DO 400 1=1»NM 

DO 305 K=1,NM 
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T+LN#=(T-I)=rf 

Тх2=7 L09 OQ 

JANILNOD 

Casu ( NV )ua*Z*(C *C)ys(OC iy 

YYA+LN#LI-I)=r 

Tx*Z=1 909 OU 

049 01 09 

3MNIJNO> 

Cual VV )a* «(C *CO)vs (OC * CON 

T+LNn#=(I-I1)=r 

92%2-:7 699 ОА 

ЗПМІ1МО> 

саж (уу) ж 2+0 67) у= (гг) 

7+1 Ме (1-71) =г 

923%2:1 899 00 

L99*999*199(DW-Z *DWW)AÀI 

SNNILNOD 

Zenl( VV De® Z4(F Pf )VÜ= (C $ r) V 

IN+«(T-I)=r 

I3*2s1 $09 OQ 

SJNNILNOD 

Cu (SB) a*Z*(CI«C*Z93* Ian Z*Z 9») V-UI«Z*Z93* |w2+23)V 
TUN®Z=1 €e9 OQ 

“09 01 09 

SNNILNOD 

С++ (98) + 2+(7 +7) у= (гг) ү 

ММ%4223 -Г 019 OQ 

295%С09%299(0М-2аОйИ04) 31 

229*119*229 (*T1-0)31 

сън (ЗЕ ) *Z2+ (PF Pivzlirer)jy 

I-I«22C 

I8N%*2=] 209 ОА 

$3903 Q3dWv15D dO 3NO SI 031У3М81 35У2 3Hl N3HM Q31V343N39 
A1SnOlA3ud XINLVW 3Hl 40 NOII1VDIJIGON V SI 9NIMO103 3Hl 
TI9*T09*T11T9(dWVv12) 31 

JNNILNOD 

ANSz=(N*1)Y 

SANS =WNS 

(4*^C)8«(^* D)8*wnSssSwns 

WN*T=f OT ОО 

°О= 05 


119 
809 


609 


819 


019 
109 


909 


999 


699 


899 


199 
609 
709 
ct? 
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31101 
3002 
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3005 
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3008 
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91 


IF (MMC#2=MC ) 3002 9300693002 
DO ' 3003 Jz1,NM 

DO 3003 L=1,K6 
I=(L-1)*NT+1 
A(I+sJ)=eə5%A(I +J) 
CONTINUE 

DO 3004 J=1,NM 

DO 3004 L=1,K6 
I=(L-1)*NT+2 
AtI»J) ze 5*A(TI»J) 
CONTINUE 

DO 3055 J=leNM 

DO 3005 K=2eNR 
A(K2+K oJ) =e5*#A(K24+KeJ) 
CONTINUE 
А(Қ5,).)):.25%А(Қ5,,)) 
СОМТІМОЕ 

GO TO 3000 

DO 3007 J=leNM 

DO 3007 і1»КІ1 
Іг((-1)%МТ41 

ACT oJ) =e5#A(T oJ) 
CONTÍNUE 

DO 3008 Jz1,NM 

DO 3008 L=1,K6 
I=(L-1)*NT+2 
A(l+J)=0e5*A(19J) 
CONTINUE 

DO 3009 Jz1»NM 

DO 3099 [=2,NR 
A(K2+2*T 9 J) =e S#A(K2+241 9 J) 
CONT INUE 
ACK5+1 oJ) =e 25*#A(K541 9J) 
CONTINUE 

DO 1000 I=1+NM 

DO 1000 J=1+NM 
B(I»J)30,0 

CONT INUE 

FINDING THE INVERSE OF MATRIX A-SQURED 
CALL MATALG (AsBeNMoNMol »DET 299) 
DO 90 Jz1,NM 

DO 90 T= 1»ММ 

A(T e«J)=06 

IF (MMC*2-MC)915,94,91 
0092 М-1»МС2»2 

DO92 Iz1»NR 

J=(N-1) *NR+2*I-1 

FA=N 

FB=MCC 

FL=FA/FB 

FNA=NNR 

FIA=I 

FIX=FIA/FNA 

FXY(J)=00 

FY(J)=020 
XFORCE(J)21,0 


GU 
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s 5 


32Ws983 

М-У3 
I*2+YN*(2-N)=f 
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0%0з(Г)АХЗ 
0%1-(Г)З3>2Ч03Х 
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YNN= YN y 
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l-IsZ*SN«&(U-N)20 
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1666 01 09 
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I+LIN#(1-1)=r 
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9x*1=1 28 00 

CC )A3«8g«CC)32803X« VV Ya * Z-s (C * n0 
ММ #1 s^ €S60Q 
ЗПМІ1МО> 
O°[l= (C )32803X 
0%0=(Г)АЗ 

0*0= (Г) АХ 

22W 4/2W 4222WNW4 
(22W) 31V013222W3 
(2W) 41V013s 2N3 
VN3/VIJsXIJ 
YNN=VNJ 

ІзУІЗ 

^ I*2xsf 
YN*T=I 666 OQ 
JANI LNOD 
O*1=(f)35404x 
0%0-(Г2А4 
0%0-(Г)АХ4 
WN4/WI4=XxI4 
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I[=VI3 

g83/V3=13 

23H=9y 

BELLE 
[*2*8N*(2-N) S£ 
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3nNI 1NO5 
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08 
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5001 
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918 
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5002 
5005 


5006 
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FL=FA/FB 

FIA=I 

FNA=NNR 

FIX=FIA/FNA 

FXY(J) 2040 

FY(J)=0+0 

XFORCE(J)=1+0 

CONTINUE 

D0954 J= 1» NM 
А(.)»/):-2.%(ААЖХҒОКСЕ(.))“88жҒУҮ(2)) 
DO 86 L*1»K1 

DO 86 K= 3eNT 

J=(L-1)*NT+K 
AtJ+J-2)=AAXXFORCE (J) 
ACJ72» J)  AA*XFORCE(CJ-2) 
IF(Q-1.)3010»3011»3010 

IF( MMC*2-MC)30125301353012 
DO 3014 L=1»X6 
I=(L-1)*NT+3 
Atl-2,1)=2.*XFORCE(I-2) 
CONTINUE 

00 3015 (-1»К6 
[=(L-1)*NT+4 
A(1-2+1)=2e*XFORCE (1-2) 
CONTINUE 
A(K5+X5+1)=2+*XFORCE(KS) 
GO TO 3010 

DO 3016 tL=1»X1 
I=(L—-1)#*NT+3 

J=(L-1)*NT+4 

At 1-2»1)=2e*XFORCE (1-2) 
A(J72»J) 29 *XFORCE(J-72) 
A(K54+19K543)=20 *#XFORCE(K5+1) 
CONTINUE 

GENERATION OF MATRIX DUE TO SHEAR 
IF (K2=NT) 500025001 95000 

IF (MCC=2) 5002 95003 95002 
GENERATION FO THE MATRIX DUE TO SH 
00918 J=4»NT»2 
ACJ723»J) 2 -FXY6 2-23) /2« 
A(J»+J-3)=-FXY(J)/2e 

00919 Ј=3, Ка», 2 
AtJ-1»J)=FXY(J-1)/2e 
AtJ»J-1)=FXY1J)/20 

GO TO 5004 
ІЕ(МС-23)5000»5005»5000 

DO 5006 J=1»NR1 
A(2*J+K7+J3)=-FXY(2*J)/20 
A(KT+J92*J)=-FXY(K7+J)/20 
A(2+2*J+K2+J)=FXY(2+2*J)/20 
A(K2+J+2+2*J)=FXY(K2+J)/20 
GO TO 5007 

LF (MMC # 2—MC ) 5008 +5009» 5008 
D0966J=NT& s KT 2 
A(J=K3»J)=-FXY(J-K3)/2-+ 

At Jo J-K3)=-FXY(J)/2e 


EAR WHEN MCC IS ODD 
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“2/(Г»2АЗАХЗ х51(1%%232У 
“2/(1ЗАХЗ s(£à *23* 00V 
Гя2%2%1М-2Ағх1| 

TyuN* T=¢88600 
“2/(ГЭАХА =(E€LN-P*F IV 
*2/(€1N-C0)AX3 =(('* € IN—[ ) V 
2ecT«x* cuz" 21 600 
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406 
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GO TO 5555 
DO1660L=2»+K1 
DO1660K=3,K4»2 
J=(L-1)*NT+K 
ACJ-K3,5 J) Z-FXYCJ-K3) /25 
AtJrJ=K3)=-=FXY(J)/20 
DO1770L=2»+K1 
DO1770K=1,NT3»2 
Jz(L-1) *NT«*K 
ACJ-NT3»J)*5 ҒХҮ(.-МТ2)/2. 
ALUJS J-NT3)7 FXY(UJM/2. 
001111 і:19КІ1 

001111 К=4,М№Т»2 
J=(L-1)*NT+K 
AtJ-39J3)=-FXY(J-3)/20 
AtJoJ=3)=-FXY(J)/20 
001114 1=1,к1 

001113 Қ-3,Қ4»2 
J=(L-1)*NT+K 
A(J-1» J) SFXYUJ-1) /2« 
ACJs J-1) 5FXYUCJ) /2« 
CONT INUE 

K33=K1-3 

IF (K33) 44413 244413 944414 
DO 24413 J=19K33 
I2J*NT 

A(1+K391)=0.0 
A(1l+1+K3)=0.0 
At1+2+KkK3-2)=0.0 
At(K3-201+2)=0.0 
CONTINUE 

IF (Q-1. )5004 2501095004 
DO 401 Lz2»K]1 
Jz(L-1)*NT*3 
AtJ-K3,5J)20, 

DO 402 L=1,K! 
J=(L-1)*NT+4 
AtJ-3,J3)=0. 

DO 403 L=1»K6 
J=(L—1)*NT+3 

ACJ-1»J) 20, 

DO 405 J=K5 ,ММ»2 
AtJ-1»J)=08 

DO 406 L=1»K6 
J=(L-1)%*NT+4 
A(J+NT39J)=008 

DO 440 J=1»NR1 
IzK2*24*2*J 

М-Қ2%2%,)-1 

А( І»М)-0. 
A(K2*25K2*3) 20, 
GENERATION OF MATRIX DUE TO FY WHEN MCC IS ODO 
IF (MMC*2-MC)2112521135»2112 
DO2114 J=2,K2»2 

AC J-1»J) :BB*FY( J-1) 
ACJS J-1) 2BB*FYCJ) 
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WN*TIszC 550% ОО 

JANI LNOD 
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9x*T=1 *00* OQ 

WN*IsC *00* OQ 
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(C*IV«S*x(C* IW 
I*1N«CI-71) 2I 

U1*1s1 €00* OQ 

ММ Тег €00% ОО 

2004 *900%* 200% ( DW—-Z # WW) SI 
000%*[00%*000%( ° T-0) 41 
0350 ЗУУ AYLSWWAS 30 
S3NI1 N3HM XIMLIVW 32803 3Hl JO NOILD3ZNYOD у SI 9NIMO1103 3H1 
(C+23)A3Aaw88#°2=([I°*[ r +Z3) V 
C#Z+LN-ZN=1 

УМ Таг 888 OQ 
2212%666%2212(%1-0)31 
(Г%2АЗАЗ3ж88:< (І1%С%2209У 
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(*¿+1N-23X= I 
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(ClAd«88si*9-C nv 
(*$Ó9-lC)A3«88s(C* 599-00v 
29%4АФСЭгГ 2<12 OQ 

2212 01 09 
(l*C)A3«88s *2Zs(£C * lO 00v 
216 A=FTETZOO 
2212%0612%2212(%1-0)31І 
(Г»АЗ3%2-(9Х-Ге%Г)У 
ММ%6АХ-Г 621200 
1212%8212%41212(%1-03341 
(ГЭАЗ%88:(%93-Г% ГУ 
(9A-P)43+88= (1 9-0 )y 
2%ММ%6іғ(921200 
*elc*Ss2l2*s2l2(0W-Z «2WW) 3I 
2212 Ol 09 
(F+21)43+99= (Za ° +Z 3) V 
(2*C)A3#88=( (+2Z3 ° Z#[ ) V 
УМ Таг Є2гтгоа 
T212*¿2l12$t1212(2W-Z+*DWW)3I 
61T12*0212*61l112(A1N-2») 3I 
(T+ AdeGGu*Zz(* Ter dy 
C*NN*$9 sC8T1T20Q 
STI2*2ZTT2*sttz(*t-O) 41 
(91-7) А-+80= (791-7) у 
(CAd3e«ggs(ooi-C* 00v 
£*NN*€»9zCOST2OQ 
(C)Adegagstt-C*o)v 
(l-C)A3«88s(C * 1-00 v 
Z2*WN*2=fP 9TTZ ОО 

STIZ OL 09 


700% 


£00% 


200% 
100% 
2212 
888 

666 

Sete 
vele 
tete 
cete 
212 
1212 
тете 
ост2 


6212 
8212 


9212 


9212 
6тт2 


EZIZ 
T2t2 
о212 
sli 
8ТІ2 
LITZ 


OSTZ 


өтт2 


етт2 


2 
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4005 
4055 


4006 


4007 
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4099 


4009 
4000 


60 
50 
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901 


1001 
3999 
1009 
11009 
11010 


11011 
1017 


1016 
8920 


8919 


A(K2*K»9J) 269 5*A(K2*K 9 J) 

CONT INUE 

A(K5eJ)=e25#A(K5e J) 

CONTINUE 

GO TO 4000 

DO 4007 J=1,NM 

DO 4007 L=1»k1 

I=(L-1)*NT+1 

ACISJ)2.4 5*A(C(I9J) 

CONTINUE 

DO 4008 Jz1,NM 

DO 4008 L31,K6 

I=(L-1)*NT+2 

ACT oJ) =e5#A(T oJ) 

CONTINUE 

DO 4009 J=1,NM 

DO 4099 I=29+NR 
AIK2+2#1,))=.5#AlK2+2#1,)) 
CONTINUE 
A(K5+19J)=.25*A(K5+19J) 
CONTINUE 

00 901 Із 1,ММ 

DO 50 K=1,NM 

SUM=O, 

DO 60 J=]1,NM 
SUMS=SUM+B( I» J)#ALJsK) 

SUM= SUMS 

Y (K)=SUM 

00320 К=1, мм 

ВІК) =Ү(К) 

CONT INUE 

DO 1001 J=1»NM 

00 1001 Із1,ММ 

A(IsJ)=B(IsJ) 

CONTINUE 

FINDING THE EIGENVALUE 

CALL EIG3(AsNMo1lsRTRoRTI 999) 
SOLVING FOR THE BUCKLING COEEFICIENT 
FNCREFLOATF (MCC) #FLOATF (NNR) *BS/(AS*#RTR(1)) 
IF(Q-14)110095,11010,11009 
GO TO 11011 

ЕМСЕ = 4, + ( ЕММЕ- 1) + (ЕМСС- 10) *BS/(AS*RTR(1)) 
IF (RUN=6«) 1017,1016,1017 
NRUN=RUN 

FORCE (NRUN) =ABSF(FNCR) 
RUN=RUN+10 

GO TO 814 

FORCE (NRUN) =ABSF(FNCR) 

PRINT 8920 

FORMAT (1H1/(40X92H 09////1///1/1/1//1111//) 
PRINT 8919+AS9»BS 

FORMAT 1(////////9+40X+14H ASPECT RATIO=E10e5»1H/»+E1005) 
DO 8883 Jz1»56 

FORCE1(Je1) =FORCE( J) 
C(J»1)=1.0 

HHX(J)zAS/XtJ) 
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C)!M/S8)a*5D2-*tasal (C) X/SV) a€9D-Can CUC)X/SV) 4232- (C)32483035 CI * C) 132303 

9*T=f ТӨӨ OQ 
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(017%023х63 Ho*XHHdLYVWHOS 6888 

$532*6888 LNIUd 

(I*$çS)I3238903=S3)D 
(0T*023s*) H**X*9*S)1vWWyO3 8886 

%32° 8888 LNIYd 

(1%%9)132804г%932 
(OLl*O23s€*9 H**X**S)1VWNOJ 16888 

£32% 888 LNIYd 

(1*€£)1323032€92 
(0OTI*023sz25 H**X** )1VWNO3 9988 

c32*98880 L1NIMd 

(1%42) І32804:222 
(0T *023s0l5» H**X**)1YWNO3 S998 

11245868 L1NIMdg 

(TI*1)13298032 U32 

(9*130*0*T1*9*132303*2)9)0V1VW YD 
3NNILNOD €888 

Zu*CUC) XHHO) &CO at CC )AHH00)s(9* 002 

Vaal (PC )AHH) E (S * 0202 

$us CCC)XHH) а (ЕГ) 

2 UUC)AHH) 8S Un * 602» 

Cue (fC) XHH) 2(2° 95D 

(7 M/SG=(7) AHH 


9801 


9800 
8950 
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13 


11 


14 


15 


16 


17 





PRINT 9801 

FORMAT (//+44X»+13H EIGENVECTORS//) 
РКІМТ 9800» ((8(.»1)».2/-1»ММ)»1-1,1) 
ҒОКМАТ( 5Е20.12) 

END 

SUBROUTINE MATALG(A»sXsNR»sNVsS IDOsDETsNACT) 
DIMENSION A(NACTsNACTO »XONACT • АСТ } 
IFCIDO) 192,1 

DO 3 I=1»NR 

DO & J=1 NR 

XCIsJ) 2060 

X(ls+1)=1.0 

NV =NR 

ОЕТ-1.0 

NR1sNR-1 

DO 5 K=1»NR1 

IR1=K+1 

РІУОТ-0.0 

DO 6 I=KeNR 

Z=ABSF(AC I 9K)) 

IF(2-PIVOT) 656,7 

PIVOT=Z 

ІРЕгі 

CONTINUE 

IF(PIVOT) 85998 

ОЕТ-0.0 

RETURN 

IF(IPR=K) 10911510 

DO 12 J=KeNR 

Z=A(IPR»J) 

A(IPRsJ)=A(K oJ) 

A(K5»J)22 

DO 13 Js1,NV 

Z=X(IPR»J) 

X(IPRsJ)=X(K +J) 

X(K»J)22 

DET=-DET ^ 
DET=DET*A(K yK) 
PIVOT=1leO/A(K 0K) 

DO 14 J=IR1 NR 
A(KeJI=A( Ke J) *PIVOT 

DO 14 I=IR1,NR 

A( T9» JJ) SACIS JJ-ACISK)*A(K9»J) 
DO 5 Jz1»NV 

IF(X(KeJ)) 1595915 
X(K»J)2X(K9 J) *PIVOT 

DO 16 I=IR1»NR 

XI» J) Z XCI» J3-ACISK)*X(K»9J) 
CONTINUE 

IFCA(NRSNR) ) 1799917 
DET=DET*#A(NR *|NR ) 
PIVOT=1e0/A(NRo NR) 

DO 18 Js1,NV 
X(NReJ)=X(NRoJ) *PI VOT 

DO 18 K=1,NR1 

I=NR=-K 
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єт 01 09 

(ON* I 18*818*AN*üN*dN**-3*T*v)839v1 IV) 
LAN*W) JON IWX=dN 

02*12*02(1-ПМ-АМ) 31 

£1 OL O9 

(QN) I18*(nN)S18*6*9 3dVI 1fdinO 3118M 
*0= (NN) 113 

(ON QN) Vs CON) HIM 

61*81*61(0N-AN) 3I 

91 01 09 
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32Vuil***9 3dV1 1ndinO 3118M 
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N*2=1 ТІ 00 

(l1*1)V=32V%91 
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6532181УИ IVI JO S301VAN39I3 ZS[IOTI #24 
(ONS TLYSYLYSW ONS VI EOI I3NLINOYGNS 
QN3 
Nun 33 
W S-=(C *I)X=(C *I)X 
(761+7) Хж(71+7 ID)vewnszwns 
ТУМ 1=7 61 OQ 
0°O=WNS 
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WRITE OUTPUT TAPE 65,6,X 
RETURN 

4 FORMAT (1H048X»+7HTRACE' =E16.8) 

5 FORMAT (11HOEIGENVALUE 12X52E20.48) 
FORMAT(1HO35Xe20HSUM OF EIGENVALUES =E16e8) 
END 

/2ALMOST TRIANGULAR (HESSENBERG) SUBROUTINE 
SUBROUTINE TRING(AsEPSeNeINTs NQ) 
DIMENSION A{(NQsNQ) > INTINO) 
WRITE OUTPUT TAPE 6»1 
Nl=N-1 
N2=N-2 
DO 21 J=leNl 
SzABSF(ACJoS J41) ) 

J1=J+1 

J22J*2 

L=J1 

NJ1zN-J1 

IF (NJ1) 155,15,6 
DO 12 K=J2eN 
T=ABSFIAULJ»K)) 
IFCT-S)12,5 12» 11 
L=K 

S=T 

CONTINUE 

IF (L-J1)13915013 
00 131 К=1,М№ 

Т=А (К, 7+1 ) 
A(KoJ+1)=A(KoL) 
A(KeLl)=T 

DO 141 K=leN 
T=A(J419K) 
A(J+1oK)=A(LoK) 
AlL»K)=T 
IF(S-EPSMINIF(ABSFIA( Jo J))»ABSFIA(J+19J+1))))16916917 
L=0 

NJ1s0 

GO TO 181 

T=A (Jo J41) 

DO 18 K#J2eN 
AtJo+K)=A(JoK)/T 
DO 20 Ічя1,М 
МаеМІМОҒ(.), 1-2) 
UzO, 

IF (NJ1) 1991907 
DO 8 KzJ2,N 
UzU*A(KoS T) *ACJsK) 
IF (M) 20,2059 

ОО 10 Кг1»М 
а0-А(К»1І)ҒА(,/ғ19К41) 
A(J+1o1)=A(J+1o]1)+U 
INT(J)sL 

INT(N) 20 

RETURN 
FORMAT(1HO48Xe22HALMOST TRIANGULAR FORM) 
END 
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YBAR=0. 

M=6 

ІҒ(ҮВАВ)27»26»27 

M=3 

DO 34 KzNUSN 

Т--А(Қ»Қ41) 

DO 34 L=loM 

SSSIGNF(1e»34 5-FLOATF(L)) 
M1sL*3*XFIXF(S) 


R-Z-XBAR*P(L,K) *YBAR* S*YP(M1,K)-FLOATF(XMODF(L-1»3) ) *P(L-15»K) 


DO 28 J=NU»K 
R=R+P(L+J)*A(K+J) 
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